# include <stdio.h>
# include <math.h>
# include <malloc.h>
# include <stdlib.h>
# include <time.h>
# define PI 3.14
# define MAX 629 //MAX是反射系数中最大的长度

/* -------------------函数定义区---------------------- */
void f(float *E); //生成雷克子波
float * h_Ip(float *pRaa1, int leng_1, float *pRaa2, int leng_2); //卷积运算 
float *f_shexi(int goat); //生成反射系数
void F_jiazao(float * pRaa, int length); //加噪
void DFT(float *x, float *xr, float *xi, int N); //x[]是信号
void f_zhengFu(float *xr, float *xi, float *Amp, float *phi, int N); //相位和振幅
void f_shexiyuanshi(void) ; //生成原始反射系数
//---------------------------大实验1----------------------------------------
void main()
{
	float E[71]; //用来存放最终雷克子波的数据
	float * fanSen = NULL; //用来指向反射系数
	float * quanJi = NULL; //指向卷积后的数据
	int i, j, K;	//用来循环的量
	FILE * fp1 = NULL; //定义一个文件类指针
	FILE * fp2 = NULL;
	FILE * fp3 = NULL;
	float Xr[699]; //DFT的的实部
	float Xi[699];	//DFT的的虚部
	float Amp[699]; //振幅谱
	float phi[699]; //相位谱
	//Xr是DFT的实部， Xi是虚部, Amp[]是增幅谱， phi[]是相位铺
	
	/* ----打开文件并判断,等下将雷克子波和反射系数的卷积值存放其中 */
	if(((fp1=fopen("quanJi","wb"))&&((fp2=fopen("Amp", "w")))&&(fp3=fopen("phi", "w"))) == NULL) 
	{
		printf("erorr!!!");
		exit(0);
	}
	/*-----------将反射系数和雷克子波的卷积写到文件----------*/
	f_shexiyuanshi();//原始数据
	f(E); //生成一个雷克子波
	for(i=0; i<200; i++)
	{	
		fanSen = f_shexi(i); //生成反射第K列的反射系数
		quanJi = h_Ip(E, 71, fanSen, MAX); //调用函数进行卷积，并将数据的地址返给quanJi
		//------------将卷积后的数据写入data文本---------------
		for(j=0; j<(MAX+71-1); j++)
		{
			fwrite(&quanJi[j], sizeof(float), 1, fp1);
		}
	}
	/*------------选择第几道卷积结果进行傅里叶变化------------*/
	printf("请输入你想让进行傅里叶变化的道数？(0~199)\n");
	scanf("%d", &K);
	fanSen = f_shexi(K);
	quanJi = h_Ip(E, 71, fanSen, MAX);
	F_jiazao(quanJi, (MAX+71-1)); //加噪
	DFT(quanJi, Xr, Xi, (MAX+71-1));
	f_zhengFu(Xr, Xi, Amp, phi, (MAX+71-1));
	for(j=0; j<(MAX+71-1); j++)
	{
		fprintf(fp2, "%f \n", Amp[j]);
		fprintf(fp3, "%f \n", phi[j]);
	}
	/*----------------关闭文件--------------*/
	if(((fclose(fp1))&&(fclose(fp2))&&(fclose(fp3))) != 0)	
	{ 
		printf("the file can not be close!");
		exit(1);
	}
	/*----------释放内存------------*/
	free(fanSen); 	
	free(quanJi);
}
/*----------雷克子波---------*/
void f(float *E) 
{
	float dt = 0.001;	//采样间隔
	int L = 71; //子波长度
	int Fm = 35;	//子波的中心频率
	int *I = NULL; //I是时间间隔点
	float *t = NULL;
	float *r = NULL;
	int i, j;
	/* 动态分配内存 */
	r = (float *)malloc(sizeof(float) * L);
	t = (float *)malloc(sizeof(float) * L);
	I = (int *)malloc(sizeof(int) * L); 
	/* 为时间间隔点赋值 从0~70 */
	for(i=0; i<L; i++)
	{
		I[i] = i;
	}
	/* 通过t = i * dt的公式计算   */
	for(i=0; i<L; i++)
	{
		t[i] = dt * I[i];
	}
	/*  通过计算r(t)的公式计算  并将值保存在R[]中 */
	for(i=0; i<L; i++)
	{
		r[i] = (exp(-1*pow(PI*Fm*t[i], 2))) * (1 - 2*pow(PI*Fm*t[i], 2));
	}
	/* --------将r[]数组变换--------------*/
	/*-----写的太傻，到时候改----*/
	for(i=35,j=0; i<=70,j<=35; i++,j++)   // 12 3 45 , 01 2 34
	{
		E[i] = r[j];
	}
	for(i=0,j=35; i<=35,j>=0; i++,j--)
	{
		E[i] = r[j];
	}
	/* 释放内存 */
	free(t);
	free(I);
	free(r);
}
/* 一维离散卷积公式 */
/* 用于计算(pRaa1[leng_1] * pRaa2[leng_2])[leng_3] */
float * h_Ip(float *pRaa1, int leng_1, float *pRaa2, int leng_2) //leng_1是1数组的长度，leng_2是2数组的长度,*pRaa3用来储存卷积后的数据
{	
	float *pRaa3 = NULL;
	int n, m; 
	int leng_3; //leng_3是存放卷积完的长度
	float k=0; //k用来临时存放数据
	pRaa3 = (float *)malloc(sizeof(float)*(leng_1+leng_2-1));//动态分配内存，！没有销毁
	leng_3 = leng_1+leng_2-1;
	/* 通过离散卷积公式计算 */
	for(n=0; n<leng_3; n++)
	{
		for(m=0; m<leng_1; m++)
		{
			if((n-m)>=0 && (n-m)<leng_2)
			{
				k += pRaa1[m] * pRaa2[n-m];
			}
		}
		pRaa3[n] = k;
		k = 0; 
	}
	return pRaa3;
}

//-----------------------------计算反射系数的函数------------------------------------------

float * f_shexi(int goat)   //goat是你选择的那一列的反射系数 一共200；   goat ======》 0 ~ 199, 而lengTh是对应反射系数数组的长度，做卷积用到 
{	
	FILE *fp = NULL;
	int i, j; //for循环中使用的循环变量
	float date[200][200]; //这是用来存放从文本读入的数据
	float t[200][3];//存放通过 2h/v = t; 来计算的时间
	float C[3] = {2000.0, 4000.0, 5000.0};//存放速度值方便计算
	int J = 0; //存放最初的j的值
	int z = 0;
	int sum = 0; //用于累加的临时变量
	int L[200]; //时间转换后每列数组的长度
	float *fShe = NULL; //反射系数 (v2 - v1)/(v2 + v1)
	//上面没有注释的大多数是单纯用来循环的量，没有别的特殊意义
	
	/* 打开文件并判断是否成功 */
	if((fp=fopen("cons","r")) == NULL)  
	{
		printf("the file can not be open");
		exit(1);
	}
	/* 通过文本的方式打开文件后，并把文本数据输入到date[][] */
	for(i=0; i<200; i++)
	{
		for(j=0; j<200; j++)
		{
			fscanf(fp, "%f", &date[j][i]);
		}
	}
	/* 下面利用公式2h/v =t 计算t 并把数据保存到t[][] 也就是计算时间，通过t*1000,找到新数组的长度*/
	/* 我把每一列速度改变一共是3段 每段的时间放在 t[][]里， 比如t[1][2] 就是第一列 第二段的时间值 */
	for(j=0; j<200; j++)
	{
		z = 0;
		J = 0;
		for(i=0; i<200; i++)
		{   
			if(date[i][j] != C[z] || i == 199) //int z = 0; C[z]是存放是速度（2000，4000，5000）方便计算
			{
				t[j][z] = (float)((i-J)*2*5) / (float)C[z]; //p是循环量 =》0~2
				z++;
				J = i; 
			}
		}
	}   
	//-------时间转换后每列数组的长度,每列分别保存到L[]里面-------

	for(i=0; i<200; i++)		

	{
		sum = 0;
		for(j=0; j<3; j++)
		{
			sum += (int)(t[i][j]*1000);
		}
		L[i] = sum;
	}
	/* 计算每一列的反射系数 */
	fShe = (float *)calloc(MAX, sizeof(float)); //L[0]-1 长度
	i = (int)(t[goat][0]*1000);
	j = (int)(t[goat][1]*1000);
	fShe[i-1] = (float)(4000-2000)/(4000+2000);
	fShe[i+j-1] = (float)(5000-4000)/(5000+4000);
	//------------关闭文件---------------

	if((fclose(fp)) != 0)

	{ 
		printf("the file can not be close!");
		exit(1);
	}
	return fShe;
}
/* ----讲数组随机赋值（加噪）------ */
void F_jiazao(float * pRaa, int leng) //length是长度
{	
    int i = 0;   
	/* 产生leng个 0.0~1.0 之间的随机浮点数 */
	/* 初始化伪随机数生成器 */
	srand((unsigned int)time(NULL));
    for(i=0; i<leng; i++)
	{	
		pRaa[i] += ((float)rand()/RAND_MAX) / (1000); //RAND_MAX是宏, 1000为倍数	
	}
}
/*--------离散傅里叶变换--------*/
void DFT(float *x, float *xr, float *xi, int N) //x[]是信号
{
	int i, j;
	for(i=0; i<N; i++)
	{
		xr[i]=0.0;
		xi[i]=0.0;	
		for(j=0; j<N; j++)
		{
			xr[i] += x[j] * cos(2*PI*i*((float)j/N));
			xi[i] -= x[j] * sin(2*PI*i*((float)j/N));
		}
	}
}
/*-----相位和振幅------*/
void f_zhengFu(float *xr, float *xi, float *Amp, float *phi, int N)
{
	int i;
	for(i=0; i<N; i++)
	{
		Amp[i] = sqrt(xr[i]*xr[i]+xi[i]*xi[i]);
		phi[i] = atan(xi[i]/xr[i]);
	}
}
void f_shexiyuanshi(void) 
{
	FILE *fp = NULL;
	int i, j; //for循环中使用的循环变量
    float date[200][200]; //这是用来存放从文本读入的数据
	FILE * fp4 = NULL;
	
	/* 打开文件并判断是否成功 */
	if((fp=fopen("cons","r"))&&(fp4=fopen("yuanshijuanji", "wb")) == NULL)  
	{
		printf("the file can not be open");
		exit(1);
	}
	/* 通过文本的方式打开文件后，并把文本数据输入到date[][] */
	for(i=0; i<200; i++)
	{
		for(j=0; j<200; j++)
		{
			fscanf(fp, "%f", &date[j][i]);
		}
	}
//===================================================
	 for(i=0; i<200; i++)
	{
		for(j=0; j<200; j++)
		{
			fwrite(&date[j][i], sizeof(float), 1, fp4);
		}
	}
	//------------关闭文件---------------

	if((fclose(fp))&&(fclose(fp4)) != 0)
	{ 
		printf("the file can not be close!");
		exit(1);
	}
}